function [x, acc] = jacobi(A,b,n,x0)

% acc = jacobi(A,x0,b,n)
%
%   Solves Ax = b using the Jacobi method.
%
%   x0 := initial x
%   n  := # of iterations

D = diag(diag(A));
Dinv = inv(D);
R = A-D;

J = Dinv * R;
c = Dinv * b;
x = x0;

acc = zeros(3, n);

for k = 1:n
    
    x = c - J*x;
    acc(1, k) = k;
    acc(2, k) = x(1);
    acc(3, k) = x(2);
    
end